############################
############################
############################
############################
if(T == T){ 
rm(list = ls() ) 
setwd("~/Dropbox/AWS Infrastructure/")
set.seed(1)
library(mvtnorm)
library(gtools)
library(colorRamps)
n_L <- 300  
n_U <- 300 
nFeat <- 2
nCat <- 2
combi_contrast <- t( combn(1:nCat, 2) )
combi_redund_alt_t <- t(combn(1:nFeat, 2))
norm1 <- nFeat * dim(combi_contrast)[1]; norm2 <- dim(combi_redund_alt_t)[1] * dim(combi_contrast)[1] 
combi_contrast_col1 <- combi_contrast[,1]; combi_contrast_col2 <- combi_contrast[,2]
combi_redund_alt_t_col1 <- combi_redund_alt_t[,1]; combi_redund_alt_t_col2 <- combi_redund_alt_t[,2]


obj_fxn_orig <- function(PrSGivenD){ 
  t_array_contrasts <- abs(PrSGivenD[,combi_contrast_col1] - PrSGivenD[,combi_contrast_col2] )
  contrasts_component <- sum( t_array_contrasts ) / norm1
  
  if(class(t_array_contrasts) == "numeric"){t_array_contrasts <- as.matrix(t_array_contrasts)}
  t_array_redun <- abs(t_array_contrasts[combi_redund_alt_t_col1,] - t_array_contrasts[combi_redund_alt_t_col2,])
  redundancy_component <- sum( t_array_redun ) / norm2
  
  return(  contrasts_component  + redundancy_component ) 
} 

#obj_fxn_alt( PrSGivenD )
#PrSGivenD <- ESGivenD_L
obj_fxn_alt <- function(PrSGivenD){ 
  t_array_contrasts <- abs(PrSGivenD[,combi_contrast_col1] - PrSGivenD[,combi_contrast_col2] )
  contrasts_component <- sum( t_array_contrasts ) / norm1
  
  if(class(t_array_contrasts) == "numeric"){t_array_contrasts <- as.matrix(t_array_contrasts)}
  t_array_redun <- abs(t_array_contrasts[combi_redund_alt_t_col1,] - t_array_contrasts[combi_redund_alt_t_col2,])
  redundancy_component <- sum( t_array_redun ) / norm2
  
  return(  redundancy_component + contrasts_component ) 
} 

contrast_ret <- function(PrSGivenD){ 
  t_array_contrasts <- abs(PrSGivenD[,combi_contrast_col1] - PrSGivenD[,combi_contrast_col2] )
  contrasts_component <- sum( t_array_contrasts ) / norm1
  
  return( contrasts_component  )
} 

redundancy_ret <- function(PrSGivenD){ 
  t_array_contrasts <- abs(PrSGivenD[,combi_contrast_col1] - PrSGivenD[,combi_contrast_col2] )
  contrasts_component <- sum( t_array_contrasts ) / norm1
  
  if(class(t_array_contrasts) == "numeric"){t_array_contrasts <- as.matrix(t_array_contrasts)}
  t_array_redun <- abs(t_array_contrasts[combi_redund_alt_t_col1,] - t_array_contrasts[combi_redund_alt_t_col2,])
  redundancy_component <- sum( t_array_redun ) / norm2
  
  return( redundancy_component  )
} 


n_rep <- 5000 
concept_drift_indicator = div_outer <- shift_outer <- redund_outer2 <- cor_outer <- rep(NA, times = n_rep)
obj_fxn_alt_outer <- obj_fxnAlt_vec <- contrasts_outer <- redund_outer <-obj_fxn_outer <- shift_vec <- mse_outer <- rep(NA, times = n_rep)
for(i in 1:n_rep){ 
  print(i)
  
  PrD_L <- round(n_L * rdirichlet(n=1, alpha=rep(2, times = nCat))[1,])
  PrD_L <- PrD_L/sum(PrD_L)
  
  PrD_U <- round(n_U * rdirichlet(n=1, alpha=rep(2, times = nCat))[1,])
  PrD_U <- PrD_U/sum(PrD_U)
  
  SigmaInput <- diag(nFeat)
  SigmaInput[] <- 0
  diag(SigmaInput) <- 1
  myDF <- nFeat+1
  #mySigma <- rWishart(n = 1, df = nFeat+1, Sigma = SigmaInput)[,,1]
  #mySigma <- 1/myDF * rWishart(n = 1, df = myDF, Sigma = SigmaInput)[,,1]
  #mySigma <- cov2cor(mySigma)
  mySigma <- diag(nFeat)
  concept_drift_indicator_i <- sample(c(0, 1), 1)
  for(asdf in 1:nCat){ 
    eval(parse(text = sprintf("Cat%s_means <- rnorm(nFeat, mean = 0, sd = sqrt(1/18))", asdf)))
    eval(parse(text = sprintf("Cat%s_means_test <- rnorm(nFeat, mean = Cat%s_means, sd = 1 * abs(Cat%s_means)*concept_drift_indicator_i)", asdf, asdf, asdf)))
  } 
  eval(parse( text =   paste("ESGivenD_L_true <- cbind(", paste(sprintf("Cat%s_means", 1:nCat), collapse = ","), ")") ) ) 
  eval(parse( text =   paste("ESGivenD_L_true_test <- cbind(", paste(sprintf("Cat%s_means_test", 1:nCat), collapse = ","), ")") ) ) 

  internalR <- 100
  inner_cor <- contrast_vec <- redundancy_vec <- div_vec <- obj_fxnAlt_vec <- obj_fxnOrig_vec <- rep(NA, times = internalR)
  estimate_mat <- matrix(NA, nrow = internalR, ncol = nCat)
  for(intR in 1:internalR){ 
  
  for(asdf in 1:nCat){ 
      eval(parse(text = sprintf("DatCat%s_L <- rmvnorm(ceiling(n_L*PrD_L[%s]), mean = Cat%s_means,sigma = mySigma)", asdf, asdf,asdf)))
      eval(parse(text = sprintf("DatCat%s_U <- rmvnorm(ceiling(n_U*PrD_U[%s]), mean = Cat%s_means_test,sigma = mySigma)", asdf, asdf,asdf)))
  } 

  eval(parse( text =   paste("BIG_TRAIN_DF <- rbind(", paste(sprintf("DatCat%s_L", 1:nCat), collapse = ","), ")") ) ) 
  eval(parse( text =   paste("BIG_TEST_DF <- rbind(", paste(sprintf("DatCat%s_U", 1:nCat), collapse = ","), ")") ) ) 
  MyMeans <- colMeans(BIG_TRAIN_DF)
  MySDs <- matrixStats::colSds(BIG_TRAIN_DF)

  #normalize 
  for(asdf in 1:nCat){ 
    eval(parse(text = sprintf("DatCat%s_L = t( (t(DatCat%s_L) - MyMeans)/MySDs)", asdf, asdf,asdf)))
  } 
  BIG_TEST_DF <- t( (t(BIG_TEST_DF) - MyMeans) / MySDs )
  
  eval(parse( text =   paste("ESGivenD_L <- cbind(", paste(sprintf("colMeans(DatCat%s_L)", 1:nCat), collapse = ","), ")") ) ) 
  ES_U <- colMeans( BIG_TEST_DF )

  my_est <- try(limSolve::lsei(A = ESGivenD_L, B=ES_U, E=matrix(nrow=1, ncol=ncol(ESGivenD_L), data=rep(1, ncol(ESGivenD_L))), 
                                F=c(1), G=diag(rep(1, ncol(ESGivenD_L))), H = rep(0, ncol(ESGivenD_L)))$X, TRUE)

  estimate_mat[intR,] <- my_est
  obj_fxnOrig_vec[intR] <- obj_fxn_orig(ESGivenD_L )
  obj_fxnAlt_vec[intR] <- obj_fxn_alt(ESGivenD_L )
  shift_vec[intR] <- mean(abs(c(ESGivenD_L_true) - c(ESGivenD_L)))
  redundancy_vec[intR] <- redundancy_ret(ESGivenD_L )
  contrast_vec[intR] <- contrast_ret(ESGivenD_L )
  } 
  
  Rez <- abs(t(estimate_mat)-PrD_U)
  ErrorVec <- colSums( Rez )
  mse_outer[i] <- mean(ErrorVec)
  div_outer[i] <- sum(abs(PrD_L-PrD_U))
  shift_outer[i] <- mean(shift_vec, na.rm = T)#concept_shift_amount#mean(abs(ESGivenD_L_true_test- ESGivenD_L_true))
  concept_drift_indicator[i] <- concept_drift_indicator_i#concept_shift_amount#mean(abs(ESGivenD_L_true_test- ESGivenD_L_true))
  obj_fxn_outer[i] <- mean(obj_fxnOrig_vec)
  obj_fxn_alt_outer[i] <- mean(obj_fxnAlt_vec)
  redund_outer[i] <- mean(redundancy_vec)
  corMat <- mySigma
  redund_outer2[i] <- mean(abs(corMat[lower.tri(corMat)]))
  contrasts_outer[i] <- mean(contrast_vec)
  cor_outer[i] <-  cor(ErrorVec, obj_fxnOrig_vec)
} 

results_mat <- data.frame(mse_outer=mse_outer, 
                     div_outer=div_outer, 
                     shift_outer=shift_outer,
                     concept_drift_indicator = concept_drift_indicator, 
                     obj_fxn_outer=obj_fxn_outer, 
                     obj_fxn_alt_outer=obj_fxn_alt_outer, 
                     redund_outer=redund_outer, 
                     redund_outer2=redund_outer2, 
                     contrasts_outer=contrasts_outer, 
                     cor_outer=cor_outer)
  write.csv(file = sprintf("./RawSimsResults_%s.csv", date() ), results_mat)
} 

results_mat <- read.csv(file = "./RawSimsResults_Tue Oct 31 19_40_10 2017.csv")

div_outer <- results_mat$div_outer
concept_drift_indicator <- results_mat$concept_drift_indicator
mse_outer <- results_mat$mse_outer
shift_outer <- results_mat$shift_outer
obj_fxn_outer <- results_mat$obj_fxn_outer
obj_fxn_alt_outer <- results_mat$obj_fxn_alt_outer
obj_fxn_outer <- results_mat$obj_fxn_outer
redund_outer <- results_mat$redund_outer
redund_outer2 <- results_mat$redund_outer2
contrasts_outer <- results_mat$contrasts_outer
cor_outer <- results_mat$cor_outer

#Fig A-1
if(T == F){ 
library(gtools)
my_q1 <- 3
x1_data <- na.omit(contrasts_outer[concept_drift_indicator==0]); x1_lab <- "Category Distinctiveness"
x1_cuts <- quantcut(x1_data, q = my_q1)
x1_levels <- tapply(x1_data, x1_cuts, function(x){median(x, na.rm = T)})

my_q2 <- 5
x2_data <- na.omit(div_outer[concept_drift_indicator==0]); x2_lab <- "Proportion Divergence"
x2_cuts <- quantcut(x2_data, q = my_q2)
x2_levels <- tapply(x2_data, x2_cuts, function(x){median(x, na.rm = T)})


x3_data <- na.omit(  mse_outer [concept_drift_indicator==0])  
mysurface <- data.frame(x = as.numeric(x1_cuts),
                        y = as.numeric(x2_cuts), 
                        z = x3_data )
mysurface <- mysurface[mysurface[,3] < quantile(mysurface[,3], probs = 0.98),]
mysurface <- mysurface[!is.na(rowSums(mysurface)),]
mysurface_rank <- mysurface

z_mat <- matrix(logical(), 
                nrow = length(unique(mysurface_rank[,1])), 
                ncol = length(unique(mysurface_rank[,2])))
my_grid <- expand.grid(c(1:my_q1), c(1:my_q2))
for(i in 1:nrow(my_grid)){
  z_mat[my_grid[i,1],my_grid[i,2]] <- mean( mysurface[mysurface[,1] == my_grid[i,1] & 
                                                mysurface[,2] == my_grid[i,2],3] ) 
}
pdf("./ContourPlot_discrim.pdf")
par(mar=c(5, 5, 4, 2) )
filled.contour(x = unique((x1_levels)),
               y = unique((x2_levels)), 
               z = z_mat, 
               nlevels=10,
               xlab = x1_lab, 
               ylab = x2_lab, 
               main = "", 
               cex.lab = 2, 
               cex.main = 2, 
               axes=TRUE,
               key.title = title(main = ("SAE")), 
               #col=hsv(h=seq(from=5/6,to=0,length=25),s=1,v=1),
               color.palette=colorRamps::ygobb)#topo.colors# heat.colors
dev.off()
} 

#Fig A-2
if(T == F){
  library(gtools)
  my_q1 <- 3
  x1_data <- na.omit(contrasts_outer[concept_drift_indicator==0]); x1_lab <- "Category Distinctiveness"
  x1_cuts <- quantcut(x1_data, q = my_q1)
  x1_levels <- tapply(x1_data, x1_cuts, function(x){median(x, na.rm = T)})
  
  my_q2 <- 3
  x2_data <- na.omit(redund_outer[concept_drift_indicator==0]); x2_lab <- "Feature Distinctiveness"
  x2_cuts <- quantcut(x2_data, q = my_q2)
  x2_levels <- tapply(x2_data, x2_cuts, function(x){median(x, na.rm = T)})
  
  x3_data <- na.omit(  mse_outer [concept_drift_indicator==0])  
  mysurface <- data.frame(x = as.numeric(x1_cuts),
                          y = as.numeric(x2_cuts), 
                          z = x3_data )
  mysurface <- mysurface[mysurface[,3] < quantile(mysurface[,3], probs = 0.98),]
  mysurface <- mysurface[!is.na(rowSums(mysurface)),]
  mysurface_rank <- mysurface
  
  z_mat <- matrix(logical(), 
                  nrow = length(unique(mysurface_rank[,1])), 
                  ncol = length(unique(mysurface_rank[,2])))
  my_grid <- expand.grid(c(1:my_q1), c(1:my_q2))
  for(i in 1:nrow(my_grid)){
    z_mat[my_grid[i,1],my_grid[i,2]] <- mean( mysurface[mysurface[,1] == my_grid[i,1] & 
                                                          mysurface[,2] == my_grid[i,2],3] ) 
  }
  pdf("./ContourPlot_joint.pdf")
  par(mar=c(5, 5, 4, 2) )
  filled.contour(x = unique((x1_levels)),
                 y = unique((x2_levels)), 
                 z = z_mat, 
                 nlevels=10,
                 xlab = x1_lab, 
                 ylab = x2_lab, 
                 main = "", 
                 cex.lab = 2, 
                 cex.main = 2, 
                 axes=TRUE,
                 key.title = title(main = ("SAE")), 
                 #col=hsv(h=seq(from=5/6,to=0,length=25),s=1,v=1),
                 color.palette=colorRamps::ygobb )#topo.colors# heat.colors
  dev.off()
  
}

#Fig B
if(T == F){ 
  library(gtools)
  my_q1 <- 4
  x1_data <- na.omit(div_outer); x1_lab <- "Proportion Divergence"
  #x1_data <- contrasts_outer; x1_lab <- "Category Distinctiveness"
  x1_cuts <- quantcut(x1_data, q = my_q1)
  x1_levels <- tapply(x1_data, x1_cuts, function(x){mean(x, na.rm = T)})
  
  my_q2 <- 3
  x2_data <- na.omit(contrasts_outer); x2_lab <- "Category Distinctiveness"
  #x2_data <- redund_outer; x2_lab <- "Feature Distinctiveness"
  x2_cuts <- quantcut(x2_data, q =my_q2)
  x2_levels <- tapply(x2_data, x2_cuts, function(x){mean(x, na.rm = T)})
  
  x3_data <- na.omit(mse_outer)
  
  x4_data <- na.omit(concept_drift_indicator)
  
  mysurface <- data.frame(x = as.numeric(x1_cuts),
                          y = as.numeric(x2_cuts), 
                          z = x3_data , 
                          alpha = x4_data )
  mysurface <- mysurface[mysurface[,3] < quantile(mysurface[,3], probs = 0.98),]
  mysurface <- mysurface[!is.na(rowSums(mysurface)),]
  
  mysurface1 <- mysurface[mysurface$alpha==0,]
  z_mat1 <- matrix(logical(), 
                  nrow = length(unique(mysurface1[,1])), 
                  ncol = length(unique(mysurface1[,2])))
  my_grid <- expand.grid(c(1:my_q1), c(1:my_q2))
  #for(i in 1:nrow(mysurface1)){z_mat[mysurface[i,1],mysurface1[i,2]] <- mysurface1[i,3]}
  for(i in 1:nrow(my_grid)){
    print(i)
    z_mat1[my_grid[i,1],my_grid[i,2]] <- mean( mysurface1[mysurface1[,1] == my_grid[i,1] & 
                                                          mysurface1[,2] == my_grid[i,2],3] ) 
  }
  
  mysurface2 <- mysurface[mysurface$alpha==1,]
  z_mat2 <- matrix(logical(), 
                   nrow = length(unique(mysurface2[,1])), 
                   ncol = length(unique(mysurface2[,2])))
  my_grid <- expand.grid(c(1:my_q1), c(1:my_q2))
  #for(i in 1:nrow(mysurface2)){z_mat[mysurface[i,1],mysurface2[i,2]] <- mysurface2[i,3]}
  for(i in 1:nrow(my_grid)){
    z_mat2[my_grid[i,1],my_grid[i,2]] <- mean( mysurface2[mysurface2[,1] == my_grid[i,1] & 
                                                            mysurface2[,2] == my_grid[i,2],3] ) 
  }
  
  if(T == T){ 
  pdf("./JointPlot.pdf", width = 10, height = 6)
  #z_mat1 <- matrix(c(0.5, 0.7,0.35, 0.55, 0.25, 0.35), ncol = 3, byrow = F)
  #z_mat2 <- matrix(c(0.70, 0.80,0.62, 0.73, 0.5, 0.65), ncol = 3, byrow = F)
  par(mfrow = c(1,2))
  par(mar=c(5, 4.5, 4, 0.7) )
  y_axis <- round(seq(0, 0.75, length.out = 5), 2)
  xbump <- -0.2
  plot(x1_levels, z_mat1[,1], ylim = c(0, max(y_axis)), type = "l", lty = 1, lwd = 3, 
       main = "Semantic Stability", cex.main = 2, yaxt = "n", col = "darkgray", 
       xlab = "Proportion Divergence", cex.lab = 1.5, ylab = "Sum of Absolute Errors (SAE)")
  axis(2, tick = T, at = y_axis, labels = y_axis)
  points(x1_levels, z_mat1[,2], col = "black", type ="l", lty = 2 , lwd = 3)
  points(x1_levels , z_mat1[,3], col = "orange", type = "l", lty = 3 , lwd = 3)
  text(par("usr")[2]-0.3*(par("usr")[2]-par("usr")[1]), max(z_mat1[nrow(z_mat1),1])+0.05, label = "Low Category Distinctiveness", col = "darkgray")
  text(par("usr")[1]+0.4*(par("usr")[2]-par("usr")[1]), z_mat1[nrow(z_mat1)/2,2]+0.10, label = "Medium Category Distinctiveness", col = "black")
  text(par("usr")[1]+0.4*(par("usr")[2]-par("usr")[1]), min(z_mat1[,ncol(z_mat1)])-0.05, label = "High Category Distinctiveness", col = "orange")
  
  par(mar=c(5, 0.7, 4, 4.5) )
  plot(x1_levels, z_mat2[,1], ylim = c(0, max(y_axis)), type = "l", lty = 1, lwd = 3, 
       main = "Semantic Change", cex.main = 2, 
       yaxt = "n",col = "darkgray", 
       xlab = "Proportion Divergence", cex.lab = 1.5, ylab = "Error")
  points(x1_levels, z_mat2[,2], col = "black", type ="l", lty = 2 , lwd = 3)
  points(x1_levels , z_mat2[,3], col = "orange", type = "l", lty = 3 , lwd = 3)
  text(par("usr")[2]-0.3*(par("usr")[2]-par("usr")[1]), max(z_mat2[nrow(z_mat2),1])+0.05, label = "Low Category Distinctiveness", col = "darkgray")
  text(par("usr")[1]+0.3*(par("usr")[2]-par("usr")[1])+0.15, z_mat2[2,2]+0.05, label = "Medium Category Distinctiveness", col = "black")
  text(par("usr")[1]+0.3*(par("usr")[2]-par("usr")[1]), min(z_mat2[,3])-0.05, label = "High Category Distinctiveness", col = "orange")
  axis(4, tick = T, at = y_axis, 
       labels = rep("", times = length(y_axis)))
  bumper_vec <- rep(0, times = length(y_axis))
  bumper_vec[nchar(y_axis) == 1] <- 0.0136
  bumper_vec[nchar(y_axis) != 1] <- 0.04
  text(par("usr")[2] + 0.06, y_axis+bumper_vec, srt=-90, adj = 0, labels = y_axis, 
       xpd = TRUE, cex = 1)
  text(par("usr")[2] + 0.14, 
       (par("usr")[4]+par("usr")[3])/2+0.40, srt=-90, adj = 0, labels = "Sum of Absolute Errors (SAE)", 
       xpd = TRUE, cex = 1.5)
  dev.off() 
  } 
  
  
}

############################